/* SPU Parallel Elliptic Curve Library - Provides Elliptic Curve Cryptographic
 * functions to the Synergistic Processing Unit target platform
 *
 * Copyright (C) 2008  Michael Annichiarico (mike.gcdeveloper@gmail.com)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
#include"field2m.h"
#include"poly_cmpeq.h"
#include"poly_cmpgt.h"

FIELD2M poly_zero = {{(ELEMENT){0ull,0ull},(ELEMENT){0ull,0ull}}};
FIELD2M poly_one = {{(ELEMENT){0ull,1ull},(ELEMENT){0ull,0ull}}};
FIELD2M poly_nist233 = {{(ELEMENT){0x400ull,1ull},(ELEMENT){0x20000000000ull,0ull}}};
/*
void poly_test() {
	printf("poly_test():\n");
	poly_print(&poly_zero);
}
*/
FIELDINLINE unsigned int poly_degree(const PFIELD2M A) {
	unsigned int total;
	vec_uint4 cnt;
	int i;
	INDEX top = poly_findtop(A);
	cnt = spu_cntlz((vec_uint4)A->v[top]);
	total = (1+(unsigned int)top)*ELEMENTBITS;
	for(i=0;i<4;++i) {
		int j=((unsigned int*)&cnt)[i];
		total-=j;
		if(j!=32) break;
	}
	if(__builtin_expect(total!=0,1)) total-=1;
	return(total);
}

FIELDINLINE unsigned int poly_degree_be(const PFIELD2M A) {
	unsigned int total;
	vec_uint4 cnt;
	int i;
	INDEX top = poly_findtop_be(A);
	cnt = spu_cntlz((vec_uint4)A->v[top]);
	total = ((unsigned int)(MAXWORDS-top))*ELEMENTBITS;
	for(i=0;i<4;++i) {
		int j=((unsigned int*)&cnt)[i];
		total-=j;
		if(j!=32) break;
	}
	if(__builtin_expect(total!=0,1)) total-=1;
	return(total);
}


FIELDINLINE INDEX poly_findtop(PFIELD2M p) {
	USELOOPVECTORS;
	ELEMENT zero = (ELEMENT)spu_splats(0);
	INDEX i=0;
//	LOOPDOWNVECTORS(p) { if(_v[_vi]==zero) {} else { i=_vi; break; } }
	LOOPDOWNVECTORSx(p, if(__builtin_expect(vector_cmpeq(_v[_vi],zero),0)) {} else {return(_vi);} )
	return(i);
}

FIELDINLINE INDEX poly_findtop_be(PFIELD2M p) {
	USELOOPVECTORS;
	ELEMENT zero = (ELEMENT)spu_splats(0);
	INDEX i=0;
//	LOOPVECTORS(p) { if(_v[_vi]==zero) {} else { i=_vi; break; } }
	LOOPVECTORSx(p, if(__builtin_expect(vector_cmpeq(_v[_vi],zero),0)) {} else {return(_vi);} )
	return(i);
}

void poly_print(PFIELD2M p) {
	USELOOPVECTORS;
	INDEX top = poly_findtop(p);
	unsigned int j;

//	LOOPVECTOR(p) { printf("%016llX %016llX",(_v[top-_vi])); if(_ei==top) break; }
	LOOPVECTORS(p) {
		for(j=0;j<NUMVECELEMS;++j) {
			VECTOROF el = ((unsigned long long)((VECTOROF*)(_v+top-_vi))[j] );
			printf("%016llx ",el);
		}
		if(_vi==top) break;
	}
}

void poly_printnice(PFIELD2M p) {
	USELOOPVECTORS;
	INDEX top = poly_findtop(p);
	int d,i;
	unsigned int j;
	int first=1;

//	printf("printnice (%p):\n",p);

	LOOPVECTORS(p) {
//		printf("  printnice vec %d (%p):\n",_vi,p);
		d=ELEMENTBITS*(top-_vi+1);
		for(j=0;j<NUMVECELEMS;++j) {
//			VECTOROF el = ((unsigned long long)((VECTOROF*)(_v+top-_vi))[j] );
		for(i=ELEMENTBITS/NUMVECELEMS-1;i>=0;--i) {
			d--;
//			printf("d: %d (top-_vi): %d j: %d i: %d\n",d,top-_vi,j,i);
//			printf("el: %0llx\n",el);
//			printf("el: %08x\n",((unsigned int*)(_v+top-_vi))[j]);
//			printf("el: %08x\n",((unsigned int*)(_v+top-_vi))[j]);
			if((unsigned long long)(((VECTOROF*)(_v+top-_vi))[j]) & (unsigned long long)(((VECTOROF)1)<<i)) {
				 if(!first) printf(" + "); printf("X^%d",d);
				first=0;
			}
		}
		}
		if(_vi==top) break;
	}
	if(first) printf("0");
}

/* From 2.41 in Rosing */
FIELDINLINE void poly_reduce_nist233(PDBLFIELD2M C) {
	USELOOPVECTORS;
	int i;
	unsigned int T;
	unsigned int* e = (unsigned int*)C->v;
	vector unsigned char endianswap = {
		0x0C,0x0D,0x0E,0x0F,
		0x08,0x09,0x0A,0x0B,
		0x04,0x05,0x06,0x07,
		0x00,0x01,0x02,0x03
	};
//	LOOPDBLVECTORSx(C, C->v[_vi]=spu_shuffle(C->v[_vi],C->v[_vi],endianswap); )
	LOOPDBLVECTORSx(C, _v[_vi]=spu_shuffle(_v[_vi],_v[_vi],endianswap); )
/*
#ifdef FIELD2M_NOUNROLL
	LOOPDBLVECTORS(C) _v[_vi]=spu_shuffle(_v[_vi],_v[_vi],endianswap);
#else
	C->v[0]=spu_shuffle(C->v[0],C->v[0],endianswap);
	C->v[1]=spu_shuffle(C->v[1],C->v[1],endianswap);
	C->v[2]=spu_shuffle(C->v[2],C->v[2],endianswap);
	C->v[3]=spu_shuffle(C->v[3],C->v[3],endianswap);
#endif
*/
	{
	for(i=15;i>=8;--i) {
		T=e[i];
		e[i-8]^=T<<23;
		e[i-7]^=T>>9;
		e[i-5]^=T<<1;
		e[i-4]^=T>>31;
	}
	T = e[7]>>9;
	e[0]^=T;
	e[2]^=T<<10;
	e[3]^=T>>22;
	e[7]&=(unsigned int)0x1FF;
	}
//	LOOPDBLVECTORSx(C, C->v[_vi]=spu_shuffle(C->v[_vi],C->v[_vi],endianswap); )
	LOOPDBLVECTORSx(C, _v[_vi]=spu_shuffle(_v[_vi],_v[_vi],endianswap); )
/*
#ifdef FIELD2M_NOUNROLL
	LOOPDBLVECTORS(C) _v[_vi]=spu_shuffle(_v[_vi],_v[_vi],endianswap);
#else
	C->v[0]=spu_shuffle(C->v[0],C->v[0],endianswap);
	C->v[1]=spu_shuffle(C->v[1],C->v[1],endianswap);
	C->v[2]=spu_shuffle(C->v[2],C->v[2],endianswap);
	C->v[3]=spu_shuffle(C->v[3],C->v[3],endianswap);
#endif
*/
}

FIELDINLINE void poly_add(PFIELD2M C, PFIELD2M A, PFIELD2M B) {
	USELOOPVECTORS;
	LOOPVECTORSx(B, C->v[_vi] = spu_xor(B->v[_vi],A->v[_vi]); )
}

FIELDINLINE void poly_shiftr1(PFIELD2M C, PFIELD2M A) {
	ELEMENT mask = makeEL{0x0000000000000000ull,1ull};
	ELEMENT carry  = makeEL{0,0};
	ELEMENT oldcarry  = makeEL{0,0};
	USELOOPVECTORS;

//	printf("shifting lower of poly: "); poly_print(A);
//	printf("\nshifting upper of poly: "); poly_print(((PFIELD2M)A)+1); printf("\n");
//#define pvec(vp) printf("vector *(" #vp ") *(%16llx) == %016llx %016llx\n",vp,((unsigned long long*)vp)[0],((unsigned long long*)vp)[1])
#define pvec(vp)

	LOOPDOWNVECTORSx(A, 
		carry = spu_and(_v[_vi],mask);
		carry=spu_rlqw(carry,-1);
		C->v[_vi] = spu_rlmaskqw(_v[_vi],-1);
		/*if(spu_cmpeq(oldcarry,mask))*/ C->v[_vi] = spu_or(C->v[_vi],oldcarry);
		oldcarry = carry;
	)
}
FIELDINLINE void poly_shiftl1(PFIELD2M C, PFIELD2M A) {
	ELEMENT mask = makeEL{0x8000000000000000ull,0};
	ELEMENT carry  = makeEL{0,0};
	ELEMENT oldcarry  = makeEL{0,0};
	USELOOPVECTORS;

//	printf("shifting lower of poly: "); poly_print(A);
//	printf("\nshifting upper of poly: "); poly_print(((PFIELD2M)A)+1); printf("\n");
//#define pvec(vp) printf("vector *(" #vp ") *(%16llx) == %016llx %016llx\n",vp,((unsigned long long*)vp)[0],((unsigned long long*)vp)[1])
#define pvec(vp)

	LOOPVECTORSx(A, 
		pvec(&_v[_vi]);
		pvec(&mask);
		carry = spu_and(_v[_vi],mask);
		pvec(&oldcarry);
		pvec(&carry);
		carry=spu_rlqw(carry,1);
		pvec(&carry);
		/*if(!spu_cmpeq(carry,makeEL{0,0})) {printf("CARRY: %016llx\n",((unsigned long long*)(&carry))[0],((unsigned long long*)(&carry))[1]);} */
		C->v[_vi] = spu_slqw(_v[_vi],1);
		/*if(spu_cmpeq(oldcarry,mask))*/ C->v[_vi] = spu_or(C->v[_vi],oldcarry);
		pvec(&C->v[_vi]);
		oldcarry = carry;
//		memcpy(&oldcarry,&carry,16);
	)
}
FIELDINLINE void poly_dblshiftl1(PDBLFIELD2M C, PDBLFIELD2M A) {
	ELEMENT mask = makeEL{0x8000000000000000ull,0};
	ELEMENT carry  = makeEL{0,0};
	ELEMENT oldcarry  = makeEL{0,0};
	USELOOPVECTORS;

//	printf("shifting lower of poly: "); poly_print(A);
//	printf("\nshifting upper of poly: "); poly_print(((PFIELD2M)A)+1); printf("\n");
//#define pvec(vp) printf("vector *(" #vp ") *(%16llx) == %016llx %016llx\n",vp,((unsigned long long*)vp)[0],((unsigned long long*)vp)[1])
#define pvec(vp)

	LOOPDBLVECTORSx(A, 
		pvec(&_v[_vi]);
		pvec(&mask);
		carry = spu_and(_v[_vi],mask);
		pvec(&oldcarry);
		pvec(&carry);
		carry=spu_rlqw(carry,1);
		pvec(&carry);
		/*if(!spu_cmpeq(carry,makeEL{0,0})) {printf("CARRY: %016llx\n",((unsigned long long*)(&carry))[0],((unsigned long long*)(&carry))[1]);} */
		C->v[_vi] = spu_slqw(_v[_vi],1);
		/*if(spu_cmpeq(oldcarry,mask))*/ C->v[_vi] = spu_or(C->v[_vi],oldcarry);
		pvec(&C->v[_vi]);
		oldcarry = carry;
//		memcpy(&oldcarry,&carry,16);
	)
}

static __inline ELEMENT _poly_makershiftmask(unsigned int count) {
	register ELEMENT tmp = (ELEMENT)spu_rlmaskqwbyte(spu_splats(~0),-((128-count)/8));
	return((ELEMENT)spu_rlmaskqw(tmp,-((128-count)&0x7)));
}
static __inline ELEMENT _poly_makelshiftmask(unsigned int count) {
	register ELEMENT tmp = (ELEMENT)spu_slqwbyte(spu_splats(~0),((128-count)/8));
	return((ELEMENT)spu_slqw(tmp,((128-count)&0x7)));
}

FIELDINLINE void poly_dblshift(PDBLFIELD2M C, PDBLFIELD2M A, signed int count) {
	if(count>0) {
	while(count>0) {
		poly_dblshiftl(C,A,count>127?127:count); A=C; count-=127;
	}
	}else{
	count*=-1;
	while(count>0) {
		poly_dblshiftr(C,A,count>127?127:count); A=C; count-=127;
	}
	}
}
FIELDINLINE void poly_shift(PFIELD2M C, PFIELD2M A, signed int count) {
	if(count>0) {
	while(count>0) {
		poly_shiftl(C,A,count>127?127:count); A=C; count-=127;
	}
	}else if(count<0){
	count*=-1;
	while(count>0) {
		poly_shiftr(C,A,count>127?127:count); A=C; count-=127;
	}
	} else if(C!=A) {
		*C=*A;
	}
}
FIELDINLINE void poly_shiftr(PFIELD2M C, PFIELD2M A, unsigned int count) {
	ELEMENT mask;
	ELEMENT carry  = makeEL{0,0};
	ELEMENT oldcarry  = makeEL{0,0};
	USELOOPVECTORS;

	mask = _poly_makershiftmask(count);

	LOOPDOWNVECTORSx(A, 
		carry = spu_and(_v[_vi],mask);
		if(__builtin_expect((128-count)/8,1)) carry=spu_slqwbyte(carry,(128-count)/8);
		if(__builtin_expect((128-count)&7,1)) carry=spu_slqw(carry,(128-count)&7);
		C->v[_vi] = spu_rlmaskqwbyte(_v[_vi],-(count/8));
		C->v[_vi] = spu_rlmaskqw(C->v[_vi],-(count&0x7));
		/*if(spu_cmpeq(oldcarry,mask))*/ C->v[_vi] = spu_or(C->v[_vi],oldcarry);
		oldcarry = carry;
	)
}
FIELDINLINE void poly_dblshiftr(PDBLFIELD2M C, PDBLFIELD2M A, unsigned int count) {
	ELEMENT mask;
	ELEMENT carry  = makeEL{0,0};
	ELEMENT oldcarry  = makeEL{0,0};
	USELOOPVECTORS;
//#define vecp(vp) printf("vector *(" #vp ") *(%16llx) == %016llx %016llx\n",vp,((unsigned long long*)vp)[0],((unsigned long long*)vp)[1])
#define vecp(vp)

	mask = _poly_makershiftmask(count);
	vecp(&mask);

	LOOPDBLDOWNVECTORSx(A,
	vecp(&oldcarry);
		carry = spu_and(_v[_vi],mask);
	vecp(&carry);
		if(__builtin_expect((128-count)/8,1)) carry=spu_slqwbyte(carry,(128-count)/8);
		if(__builtin_expect((128-count)&7,1)) carry=spu_slqw(carry,(128-count)&7);
	vecp(&carry);
		C->v[_vi] = spu_rlmaskqwbyte(_v[_vi],-(count/8));
		C->v[_vi] = spu_rlmaskqw(C->v[_vi],-(count&0x7));
	vecp(&C->v[_vi]);
		/*if(spu_cmpeq(oldcarry,mask))*/ C->v[_vi] = spu_or(C->v[_vi],oldcarry);
		oldcarry = carry;
	)
}
FIELDINLINE void poly_shiftl(PFIELD2M C, PFIELD2M A, unsigned int count) {
	ELEMENT mask;
	ELEMENT carry  = makeEL{0,0};
	ELEMENT oldcarry  = makeEL{0,0};
	USELOOPVECTORS;

	mask = _poly_makelshiftmask(count);
//	printf("count = %d : poly_shiftl mask= %016llx %016llx\n",count, ((VECTOROF*)&mask)[0] ,((VECTOROF*)&mask)[1]) ;

	LOOPVECTORSx(A,
		carry = spu_and(_v[_vi],mask);
//	printf("poly_shiftl carry= %016llx %016llx\n", ((VECTOROF*)&carry)[0] ,((VECTOROF*)&carry)[1]) ;

		if(__builtin_expect((128-count)/8,1)) carry=spu_rlmaskqwbyte(carry,-((128-count)/8));
//	printf("poly_shiftl carry= %016llx %016llx\n", ((VECTOROF*)&carry)[0] ,((VECTOROF*)&carry)[1]) ;
		if(__builtin_expect((128-count)&7,1)) carry=spu_rlmaskqw(carry,-((128-count)&7));
//	printf("poly_shiftl carry= %016llx %016llx\n", ((VECTOROF*)&carry)[0] ,((VECTOROF*)&carry)[1]) ;
//	printf("poly_shiftl C= %016llx %016llx\n", ((VECTOROF*)(C->v+_vi))[0],((VECTOROF*)(C->v+_vi))[1]);
		C->v[_vi] = spu_slqwbyte(_v[_vi],(count/8));
//	printf("poly_shiftl C= %016llx %016llx\n", ((VECTOROF*)(C->v+_vi))[0],((VECTOROF*)(C->v+_vi))[1]);
		C->v[_vi] = spu_slqw(C->v[_vi],(count&0x7));
//	printf("poly_shiftl C= %016llx %016llx\n", ((VECTOROF*)(C->v+_vi))[0],((VECTOROF*)(C->v+_vi))[1]);

		/*if(spu_cmpeq(oldcarry,mask))*/ C->v[_vi] = spu_or(C->v[_vi],oldcarry);
//	printf("poly_shiftl oldcarry= %016llx %016llx\n", ((VECTOROF*)&oldcarry)[0] ,((VECTOROF*)&oldcarry)[1]) ;
//	printf("poly_shiftl C= %016llx %016llx\n", ((VECTOROF*)(C->v+_vi))[0],((VECTOROF*)(C->v+_vi))[1]);
		oldcarry = carry;
	)
}
FIELDINLINE void poly_dblshiftl(PDBLFIELD2M C, PDBLFIELD2M A, unsigned int count) {
	ELEMENT mask = makeEL{0x8000000000000000ull,0};
	ELEMENT carry  = makeEL{0,0};
	ELEMENT oldcarry  = makeEL{0,0};
	USELOOPVECTORS;

	mask = _poly_makelshiftmask(count);

	LOOPDBLVECTORSx(A,
		carry = spu_and(_v[_vi],mask);

		if(__builtin_expect((128-count)/8,1)) carry=spu_rlmaskqwbyte(carry,-((128-count)/8));
		if(__builtin_expect((128-count)&7,1)) carry=spu_rlmaskqw(carry,-((128-count)&7));
		C->v[_vi] = spu_slqwbyte(_v[_vi],(count/8));
		C->v[_vi] = spu_slqw(C->v[_vi],(count&0x7));

		/*if(spu_cmpeq(oldcarry,mask))*/ C->v[_vi] = spu_or(C->v[_vi],oldcarry);
		oldcarry = carry;
	)
}

#define POLY_MULT_W 4
#define POLY_MULT_L 4
void poly_multw(PDBLFIELD2M C, unsigned char poly8, PFIELD2M B) {
	int k;
	unsigned char j;
	ELEMENT ke = (ELEMENT){0x0ull,0x8ull};
	USELOOPVECTORS;
	LOOPDBLVECTORSx(C, _v[_vi]=(ELEMENT)spu_splats(0); )

	for(k=POLY_MULT_W-1;k>=0;--k) {
		j=poly8&(1<<k);
			if(j) {//if(vector_cmpeq(spu_and(A->v[_vi],ke),ke)) {
//		for(_vi=0;_vi<MAXWORDS;++_vi) {
				poly_add((PFIELD2M)&(C->v[0]),(PFIELD2M)&(C->v[0]),B);
//		}
			}
		if(__builtin_expect(k,1)) {/* shift left C */
			poly_dblshiftl1(C,C);
		}
	}
//	printf("poly_multw yields: "); poly_print(C); printf(" :: "); poly_print(((PFIELD2M)C)+1); printf("\n");
}
/* using 2.36 */
FIELDINLINE void poly_mult(PDBLFIELD2M C, PFIELD2M A, PFIELD2M B) {
	int k;
	unsigned int j;
	unsigned char poly8;
//	DBLFIELD2M Bvi[1<<POLY_MULT_W/POLY_MULT_L][POLY_MULT_L];
	DBLFIELD2M Bu[1<<POLY_MULT_W];
	PDBLFIELD2M Bup;
	ELEMENT ke = (ELEMENT){0x8000000000000000ull,0};
	USELOOPVECTORS;
	LOOPDBLVECTORSx(C, _v[_vi]=(ELEMENT)spu_splats(0); )

	//Bu[0]=C[0];
	LOOPDBLVECTORSx(Bu, _v[_vi]=(ELEMENT)spu_splats(0); )
	//Bu[1]=B;
	LOOPVECTORSx((Bu+1), _v[_vi]=B->v[_vi]; )
	for(k=2;k<(1<<POLY_MULT_W);++k) { poly_multw(&Bu[k],(unsigned char)k,B);/* poly_reduce_nist233(&Bu[k]);*/ }

/*	printf("new mult test\n"); {
		int ii;
		printf("A = "); poly_printnice(A);
		printf("\nB = "); poly_printnice(B); printf("\n");
		for(ii=0;ii<(1<<POLY_MULT_W);++ii)
		printf("Bu[%d]: ",ii);
		poly_print((PFIELD2M)(Bu+ii));printf(" :: ");
		poly_print(((PFIELD2M)(Bu+ii))+1);printf("\n");
	}*/

	for(k=(ELEMENTBITS/(POLY_MULT_W))-1;k>=0;--k) {
		for(_vi=0;_vi<MAXWORDS;++_vi) {
			/*if(vector_cmpeq(spu_and(A->v[_vi],ke),ke)) {
				poly_add((PFIELD2M)&(C->v[_vi]),(PFIELD2M)&(C->v[_vi]),B);
			}*/
			//Bup = &Bu[spu_extract((vec_uchar16)(B->v[k/32]),15-(k%16))];
#if POLY_MULT_W == 4
			poly8 = spu_extract((vec_uchar16)(A->v[_vi]),15-(POLY_MULT_W*k/8));
			poly8>>=4*(k&1);
			poly8&=0xF;
#else
#error "POLY_MULT_W MUST EQUAL 4"
			poly8 = spu_extract((vec_ushort8)(A->v[_vi]),7-(POLY_MULT_W*k/16));
			{static int shiftamount[]= {0,6,4,2};
			poly8>>=shiftamount[((POLY_MULT_W*k)&7)];
			poly8&=0x3F;}
#endif
			//poly_add((PFIELD2M)&(C->v[_vi]),(PFIELD2M)&(C->v[_vi]),&Bu[spu_extract((vec_uchar16)(B->v[k/16]),15-(k%16))]);
			poly_add((PFIELD2M)&(C->v[_vi]),(PFIELD2M)&(C->v[_vi]),&Bu[poly8]);
			
		}
		if(__builtin_expect(k,1)) {/* shift left C */
			poly_dblshiftl(C,C,POLY_MULT_W);
		}
	}
}

/* using 2.35 */
FIELDINLINE void poly_mult_old(PDBLFIELD2M C, PFIELD2M A, PFIELD2M B) {
	int k;
	unsigned int j;
	ELEMENT ke = (ELEMENT){0x8000000000000000ull,0};
	USELOOPVECTORS;
	LOOPDBLVECTORSx(C, _v[_vi]=(ELEMENT)spu_splats(0); )

	for(k=ELEMENTBITS-1;k>=0;--k) {
//		LOOPVECTORSx(A, //
		for(_vi=0;_vi<MAXWORDS;++_vi) {
			//if kth bit of A[j] then add B to C{j}
//			if(A->e[(j*NUMVECELEMS)+(k/sizeof(VECTOROF))] & (1<<(k%sizeof(VECTOROF)))) {
			
//#define pvec(vp) printf("vector *(" #vp ") *(%16llx) == %016llx %016llx\n",vp,((unsigned long long*)vp)[0],((unsigned long long*)vp)[1])
#define pvec(vp)
			pvec(&A->v[_vi]);
			pvec(&ke);
			if(vector_cmpeq(spu_and(A->v[_vi],ke),ke)) {
//				printf("adding!\n");
				poly_add((PFIELD2M)&(C->v[_vi]),(PFIELD2M)&(C->v[_vi]),B);
			}
		}//)
		if(__builtin_expect(k,1)) {/* shift left C */
			poly_dblshiftl1(C,C);
		}
		ke = spu_rlmaskqw(ke,-1); /* shift right 1 bit */
	}
}

FIELDINLINE void poly_mod_mult_dbl(PDBLFIELD2M D, PFIELD2M A, PFIELD2M B) {
	poly_mult(D,A,B);
	poly_reduce_nist233(D);
}

FIELDINLINE void poly_mod_mult(PFIELD2M C, PFIELD2M A, PFIELD2M B) {
	DBLFIELD2M D;
	poly_mod_mult_dbl(&D,A,B);
	*C=*(PFIELD2M)&D;
}

/* 2.33 */
/*
void poly_mod_mult(PFIELD2M C, PFIELD2M A, PFIELD2M B) {
	int i;
	ELEMENT ke = (ELEMENT){0x0ull,1};
	DBLFIELD2M b;
	ELEMENT mask = (ELEMENT){0ull,1ull};
	*(PFIELD2M)&b = *B;
	*(((PFIELD2M)&b)+1) = poly_zero;
	if(spu_and(A->v[0],mask)==mask) *C=*B; else *C=poly_zero;

	for(i=1;i<=127;++i) {
		ke = spu_slqw(ke,1);
//		printf("b= "); poly_printnice((PFIELD2M)&b);printf("\n");
		poly_dblshiftl1(&b,&b);
		poly_reduce_nist233(&b);
//		printf("b= bz mod f "); poly_printnice((PFIELD2M)&b);printf("\n");
		if(spu_and(A->v[0],ke)==ke) poly_add(C,C,(PFIELD2M)&b);
	}
	ke = (ELEMENT){0x0ull,1};
	for(i=128;i<=233;++i) {
//		printf("b= "); poly_printnice((PFIELD2M)&b);printf("\n");
		poly_dblshiftl1(&b,&b);
		poly_reduce_nist233(&b);
//		printf("b= bz mod f "); poly_printnice((PFIELD2M)&b);printf("\n");
		if(spu_and(A->v[1],ke)==ke) poly_add(C,C,(PFIELD2M)&b);
		ke = spu_slqw(ke,1);
	}
}
*/

FIELDINLINE void poly_sqr(PDBLFIELD2M D, PFIELD2M A) {
	USELOOPVECTORS;
	vector unsigned char expandtable = (vec_uchar16){
		0x00,0x01,0x04,0x05,
		0x10,0x11,0x14,0x15,
		0x40,0x41,0x44,0x45,
		0x50,0x51,0x54,0x55
	};
	vector unsigned char alternateselect = {
		0x00,0x10,
		0x01,0x11,
		0x02,0x12,
		0x03,0x13,
		0x04,0x14,
		0x05,0x15,
		0x06,0x16,
		0x07,0x17
	};
	vector unsigned char mask1 = (vec_uchar16)spu_splats(0x0F0F0F0F);
	vector unsigned char mask2 = (vec_uchar16)spu_splats(0xF0F0F0F0);
	vec_uchar16 e1,e2;
//	vec_uchar16 res;

/*
#define pvec(vp) printf("vector *(" #vp ") *(?) == %016llx %016llx\n",((unsigned long long*)vp)[0],((unsigned long long*)vp)[1])
#define dspu_shuffle(a,b,c) spu_shuffle(a,b,c); res=(vec_uchar16)spu_shuffle(a,b,c);  \
	{vec_ullong2 A=(vec_ullong2)a, B= (vec_ullong2)b,C=(vec_ullong2)c;  pvec(a,&A); pvec(b,&B); pvec(c,&C); pvec(shuffle,&res); } 
*/
		
	LOOPDOWNVECTORSx(A,
		e2 = spu_shuffle(
					(vec_uchar16)expandtable,
					(vec_uchar16)mask2,
					(vec_uchar16) spu_and(_v[_vi],(vec_ullong2)mask1));
		e1 = spu_shuffle(
					(vec_uchar16)expandtable,
					(vec_uchar16)mask1,
					(vec_uchar16)spu_rlmaskqw(spu_and(_v[_vi],(vec_ullong2)mask2),-4));
		D->v[_vi*2+1]=(vec_ullong2)spu_shuffle(e1,e2,alternateselect);
		D->v[_vi*2+0]=(vec_ullong2)spu_shuffle(e1,e2,(vec_uchar16)spu_add((vec_ushort8)alternateselect,(vec_ushort8)spu_splats(0x08080808)));
	)
}

FIELDINLINE void poly_mod_sqr_dbl(PDBLFIELD2M D, PFIELD2M A) {
	poly_sqr(D,A);
	poly_reduce_nist233(D);
}

FIELDINLINE void poly_mod_sqr(PFIELD2M C, PFIELD2M A) {
	DBLFIELD2M D;
	poly_mod_sqr_dbl(&D, A);
	*C=*(PFIELD2M)&D;
}

/* 2.49 */
FIELDINLINE void poly_inv_mod_2dot49(PFIELD2M U, PFIELD2M A) {
	ELEMENT mask = (ELEMENT){0x0ull,1ull};
	FIELD2M V = poly_nist233;
	FIELD2M G1=poly_one;
	FIELD2M G2=poly_zero;
	if(U!=A) *U=*A;
/*
#define printpoly(p) \
	printf("Poly " #p ": "); poly_print(p); printf("\n");
*/

//	while(!poly_cmpeq(U,&poly_one) && !poly_cmpeq(&V,&poly_one)) {
	while(1) {
		if(__builtin_expect(poly_cmpeq(U,&poly_one),0)) {
//			printf("U==1\n");
//			printpoly(&G1);
			*U=G1;
			return;
		} else if(__builtin_expect(poly_cmpeq(&V,&poly_one),0)) {
//			printf("V==1\n");
//			printpoly(&G2);
			*U=G2;
			return;
		}
#if 0
		printf("top: k==%d\n",k);
		vecp(U);
		vecp(&V);
		vecp(&G1);
		vecp(&G2);
		printf("U loop\n");
#endif
		while(1) {
		ELEMENT tmp=spu_and(U->v[0],mask);
		vecp(&tmp);
		if(vector_cmpeq(tmp,mask)) break;
			//shift left u
//		printpoly(U);
			poly_shiftr1(U,U);
//		printpoly(U);
//		printpoly(&G1);
			// if z | g1 then g1 >> 1, else g1 = (g1+f)>>1
			tmp = spu_and(G1.v[0],mask);
			if(vector_cmpeq(tmp,mask))  { poly_add(&G1,&G1,&poly_nist233); poly_shiftr1(&G1,&G1); }
			else {  poly_shiftr1(&G1,&G1); }
//		printpoly(&G1);
		}
//		printf("V loop\n");
		while(1) {
		ELEMENT tmp=spu_and(V.v[0],mask);
		vecp(&tmp);
		if(vector_cmpeq(tmp,mask)) break;
			//shift left v, 
//		printpoly(&V);
			poly_shiftr1(&V,&V);
//		printpoly(&V);
//		printpoly(&G2);
			tmp = spu_and(G1.v[0],mask);
			if(vector_cmpeq(tmp,mask))  { poly_add(&G2,&G2,&poly_nist233); poly_shiftr1(&G2,&G2); }

			else { poly_shiftr1(&G2,&G2); }
//		printpoly(&G2);
		}
		//if(poly_degree(U) > poly_degree(&V)) {
		if(poly_cmpgt(U,&V)) {
//			printf("U>V\n");
			poly_add(U,U,&V);
			poly_add(&G1,&G1,&G2);
		}else{
//			printf("U <= V\n");
			poly_add(&V,U,&V);
			poly_add(&G2,&G1,&G2);
		}

	}
}


/* 2.50 */
/*
void poly_inv_mod(PFIELD2M U, PFIELD2M A) {
	USELOOPVECTORS;
	ELEMENT mask = (ELEMENT){0x0ull,1ull};
	FIELD2M V = poly_nist233;
	DBLFIELD2M G1;
	DBLFIELD2M G2;
	INDEX k=0;
	*((PFIELD2M)&G1)=poly_one;
	*(((PFIELD2M)&G1)+1)=poly_zero;
	*((PFIELD2M)&G2)=poly_zero;
	*(((PFIELD2M)&G2)+1)=poly_zero;
	if(U!=A) *U=*A;
#define printpoly(p) \
	printf("Poly " #p ": "); poly_print(p); printf("\n");

	while(1) {
		printf("top: k==%d\n",k);
		vecp(U);
		vecp(&V);
		vecp(&G1);
		vecp(&G2);
		printf("U loop\n");
		while(1) {
		ELEMENT tmp=spu_and(U->v[0],mask);
		vecp(&tmp);
		if(tmp==mask) break;
			//shift right u, shift left g2, k++
		printpoly(U);
			poly_shiftr1(U,U);
		printpoly(U);
		printpoly(&G2);
			poly_dblshiftl1(&G2,&G2);
		printpoly(&G2);
			++k;
		printf("k = %d\n",k);
		}
		printf("V loop\n");
		while(1) {
		ELEMENT tmp=spu_and(V.v[0],mask);
		vecp(&tmp);
		if(tmp==mask) break;
			//shift right v, shift left g1, k++
		printpoly(&V);
			poly_shiftr1(&V,&V);
		printpoly(&V);
		printpoly(&G1);
			poly_dblshiftl1(&G1,&G1);
		printpoly(&G1);
			++k;
		}
		if(poly_degree(U) > poly_degree(&V)) {
			printf("U>V\n");
			poly_add(U,U,&V);
			poly_add(((PFIELD2M)&G1)+0,((PFIELD2M)&G1)+0,((PFIELD2M)&G2)+0);
			poly_add(((PFIELD2M)&G1)+1,((PFIELD2M)&G1)+1,((PFIELD2M)&G2)+1);
		}else{
			printf("U <= V\n");
			poly_add(&V,U,&V);
			poly_add(((PFIELD2M)&G2)+0,((PFIELD2M)&G1)+0,((PFIELD2M)&G2)+0);
			poly_add(((PFIELD2M)&G2)+1,((PFIELD2M)&G1)+1,((PFIELD2M)&G2)+1);
		}

		if(poly_cmpeq(U,&poly_one)) {
			printf("U==1\n");
			printpoly(&G1);
			poly_dblshiftr(&G1,&G1,k);
			poly_reduce_nist233(&G1);
			*U=*(PFIELD2M)&G1;
			printpoly(U);
			return;
		} else if(poly_cmpeq(&V,&poly_one)) {
			printf("V==1\n");
			printpoly(&G2);
			poly_dblshiftr(&G2,&G2,k);
			poly_reduce_nist233(&G2);
			*U=*(PFIELD2M)&G2;
			printpoly(U);
			return;
		}
	}
}
*/

/* 2.48 */
FIELDINLINE void poly_inv_mod_2dot48(PFIELD2M U, PFIELD2M A) {
	FIELD2M V = poly_nist233;
	FIELD2M G1;
	FIELD2M G2;
	FIELD2M T;
	PFIELD2M u=U,v=&V,g1=&G1,g2=&G2;
	*g1=poly_one;
	*g2=poly_zero;
	int j;
	if(U!=A) *U=*A;
/*
#define printpoly(p) \
	printf("Poly " #p ": "); poly_print(p); printf("\n");
*/
#define printpoly(p)

	while(1) {
		if(__builtin_expect(poly_cmpeq(u,&poly_one),0)) {
//			printf("U==1\n");
			printpoly(g1);
			*U=*(PFIELD2M)g1;
			return;
		}
		vecp(u);
		vecp(v);
		vecp(g1);
		vecp(g2);
//		printf("deg(u) == %d\n",poly_degree(u));
//		printf("deg(v) == %d\n",poly_degree(v));
		j=poly_degree(u)-poly_degree(v);
//		printf("j=%d\n",j);
		if(j<0) {
			PFIELD2M t;
			t=v; v=u; u=t;
			t=g1; g1=g2; g2=t;
			j=-1*j;
		}
		poly_shift(&T,v,j);
		poly_add(u,u,&T);
		poly_shift(&T,g2,j);
		poly_add(g1,g1,&T);

	}
}


__inline void poly_inv_mod(PFIELD2M U, PFIELD2M A) {
	poly_inv_mod_2dot48(U,A);
//	poly_inv_mod_2dot49(U,A);
}


